Analisis de Datos Omicos en COVID
Alex Sanchez
Enero 2022
1 Introducción
Este estudio presenta un análisis básico de RNA-Seq con énfasis en las posibilidades de visualización de datos y resultados^.
En su versión actual los gráficos son correctos pero poco elaborados. Para mejorarlos se utilizaran gráficos con ggplot2 utilizando como modelo algunas de las propuestas del documento: https://github.com/UofABioinformaticsHub/DataVisualisaton_BIS2016.git
1.1 Infraestructura informática para el análisis
Este estudio se llevará a cabo usando R/Bioconductor por lo que es preciso tener instaladas un conjunto de librerías. Esto puede hacerse siguiendo el procedimiento descrito a continuación.
El código que se presenta debería ejecutarse tan sólo una vez!
if (!require(BiocManager)) {
install.packages("BiocManager", dep = TRUE)
}
installifnot <- function(pckgName, BioC = TRUE) {
if (BioC) {
if (!require(pckgName, character.only = TRUE)) {
BiocManager::install(pckgName)
}
} else {
if (!require(pckgName, character.only = TRUE)) {
install.packages(pckgName, dep = TRUE)
}
}
}
installifnot("limma")
installifnot("edgeR")
installifnot("org.Hs.eg.db")
installifnot("clusterProfiler")
installifnot("dplyr", BioC = FALSE)
installifnot("gplots", BioC = FALSE)
installifnot("ggvenn", BioC = FALSE)
installifnot("pheatmap")
installifnot("prettydoc", BioC = FALSE)
installifnot("ggnewscale", BioC = FALSE)
# Pels grafics que empren ggplot2
installifnot("locfit", BioC = FALSE)
installifnot("magrittr", BioC = FALSE)
installifnot("statmod", BioC = FALSE)
# Pel HeatMap Interactiu
installifnot("heatmaply")
# Pel VolcanoPlot interactiu
installifnot("here")
installifnot("janitor") # Cleaning column names
installifnot("scales") # Transform axis scales
installifnot("ggrepel")Puede resultar útil, aunque no es para nada imprescindible, disponer de, al menos dos, directorios específicos: - datos (o un nombre similar) donde guardar y de donde cargar los datos. - resultso (o un nombre similar) donde escribir los resultados.
if (!dir.exists("datos")) dir.create("datos")
if (!dir.exists("results")) dir.create("results")
if (!dir.exists("figures")) dir.create("figures")2 Los datos para el análisis
Este análisis se basa en un estudio publicado recientemente (Arunachalam 2020) que investigaba la respuesta inmune a la infección con SARS-COV-2 desde un a perspectiva de biología de sistemas utilizando tecnología de secuenciación, RNA-Seq.
Los datos se han depositado en el repositorio “Gene Expression Omnibus (GEO)” con el identificador GSE152418.
2.1 Lectura de la matriz de contajes
Los datos de contajes se han descargado desde l repositorio “GEO” al archivo Rawcounts.csv que se utilizará como punto de partida para el análisis.
counts <- read.csv("datos/RawCounts.csv", row.names = 1)
selectedCounts <- as.matrix(counts)En la misma web se dispone de información los grupos y otras covariables o variables auxiliares. Con éstos se crea un objeto (habitualmente denominado “targets”) que debe estar sincronizado con el anterior, es decir, cuyas filas deben corresponderse con las columnas de la matriz de datos.
sampleNames <- colnames(selectedCounts)
grupos <- c(rep("COVID", 17), rep("SANO", 17))
colores = c(rep("red", 17), rep("blue", 17))
selectedTargets <- data.frame(samples = sampleNames, group = grupos, cols = colores)
rownames(selectedTargets) <- selectedTargets[, 1]3 Preprocesado de los datos
3.1 Estandarización de los contajes
Los datos de secuenciación pueden estar “desbalanceados” en el sentido que distintas muestras pueden contener un número distinto de secuencias, lo que puede inducir erróneamente a pensar que un gen se expresa más en una muestra que en otra, cuando esto se deba a esta diferencia global.
Esto puede evitarse expresando los contajes como “CPMs” es decir “counts per million,” lo que no modificará los resultados de comparaciones posteriores, pero hará que las muestras sean comparables en número , lo que es útil y necesario para los análisis posteriores.
library(edgeR)
selectedCounts[1:5, 1:6] COV145 COV147 COV149 COV150 COV151 COV152
ENSG00000223972 0 0 0 0 0 0
ENSG00000227232 1 0 3 16 1 2
ENSG00000278267 3 1 2 0 3 0
ENSG00000243485 2 0 0 1 0 0
ENSG00000284332 0 0 0 0 0 0apply(selectedCounts, 2, sum) COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154
15048983 13946337 13520020 18863644 13659353 13279390 12967503 14300187
COV155 COV156 COV157 COV175 COV176 COV177 COV178 COV179
17192277 15626012 17041245 13415720 13679816 16092084 17379155 13672179
COV180 HEA061 HEA062 HEA063 HEA064 HEA065 HEA066 HEA067
13722892 17003281 17359694 15570080 16172199 19420885 18264376 16398893
HEA068 HEA069 HEA070 HEA071 HEA181 HEA182 HEA183 HEA184
14895860 14939005 20458479 17820834 14973734 14280768 18297901 17177433
HEA185 HEA186
16467873 16860317 counts.CPM <- cpm(selectedCounts)
counts.CPM[1:5, 1:6] COV145 COV147 COV149 COV150 COV151 COV152
ENSG00000223972 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991 0.1506093
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973 0.0000000
ENSG00000243485 0.13289935 0.00000000 0.0000000 0.05301203 0.00000000 0.0000000
ENSG00000284332 0.00000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000apply(counts.CPM, 2, sum)COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155 COV156 COV157
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
COV175 COV176 COV177 COV178 COV179 COV180 HEA061 HEA062 HEA063 HEA064 HEA065
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
HEA066 HEA067 HEA068 HEA069 HEA070 HEA071 HEA181 HEA182 HEA183 HEA184 HEA185
1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06 1e+06
HEA186
1e+06 Una vez los datos estan como CPMs, se procede a filtrarlos,
3.2 Filtraje de genes poco expresados
Los genes con recuentos muy bajos en todas las librerías (es decir en todas las muestras) proporcionan poca evidencia de expresión diferencial, por lo que es habitual eliminar aquellos genes que, o bien varían muy poco entre grupos, o bien presentan poca o nula expresión en la mayoría de las muestras.
En este caso, siguiendo un criterio habitual, se opta por conservar únicamente aquellos genes que presentan algún valor en, al menos, tres muestras de cada grupo.
thresh <- counts.CPM > 0
keep <- (rowSums(thresh[, 1:10]) >= 3) & (rowSums(thresh[, 11:20]) >= 3)
counts.keep <- counts.CPM[keep, ]
dim(counts.CPM)[1] 60683 34dim(counts.keep)[1] 33039 34Aunque no sea más que un ejemplo basta con ver los dos primeros genes para comprobar como el primero no cumple la condición, en el grupo “SANO,” mientras que los siguientes sí que la cumplen. Por lo tanto, al filtrar desaparece el primer gen de la matriz filtrada, pero no los dos siguientes.
round(head(counts.CPM), 3) COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155
ENSG00000223972 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.00 0.116
ENSG00000227232 0.066 0.000 0.222 0.848 0.073 0.151 1.388 0.42 1.163
COV156 COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061
ENSG00000223972 0.000 0.117 0.000 0.000 0.000 0.058 0.000 0.000 0.000
ENSG00000227232 0.000 0.117 0.000 0.950 1.989 0.173 0.146 0.656 0.353
HEA062 HEA063 HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070
ENSG00000223972 0.000 0.000 0.062 0.000 0.000 0.000 0.000 0.000 0.000
ENSG00000227232 0.000 0.000 0.062 0.103 0.055 0.000 0.067 0.067 0.049
HEA071 HEA181 HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000223972 0.000 0.0 0.00 0.000 0.000 0.000 0.000
ENSG00000227232 0.112 0.2 0.28 0.055 0.000 0.000 0.415
[ reached getOption("max.print") -- omitted 4 rows ]round(head(counts.keep), 3) COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155
ENSG00000227232 0.066 0.000 0.222 0.848 0.073 0.151 1.388 0.420 1.163
ENSG00000278267 0.199 0.072 0.148 0.000 0.220 0.000 0.000 0.070 0.174
COV156 COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061
ENSG00000227232 0.000 0.117 0.000 0.950 1.989 0.173 0.146 0.656 0.353
ENSG00000278267 0.128 0.117 0.224 0.219 0.435 0.000 0.219 0.000 0.235
HEA062 HEA063 HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070
ENSG00000227232 0.000 0.000 0.062 0.103 0.055 0.000 0.067 0.067 0.049
ENSG00000278267 0.288 0.193 0.495 0.154 0.493 0.427 0.134 0.469 0.098
HEA071 HEA181 HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000227232 0.112 0.200 0.28 0.055 0.000 0.000 0.415
ENSG00000278267 0.337 0.200 0.14 0.328 0.699 0.304 0.178
[ reached getOption("max.print") -- omitted 4 rows ]3.3 Uso de clases específicas para manejar los datos
Cuando se trabaja con distintos objetos referidos a unos mismos datos, como la matriz de contajes y el objeto “targets,” es útil disponer de contenedores que permitan trabajar con todos ellos a la vez, lo que no sólo facilita el trabajo sino que ayuda a evitar “desincronizaciones.”
Éste es el caso de la clase ExpressionSet habitualmente utilizada con microarrays o de la clase que la generaliza, llamada SummarizedExperiment.
Para datos de contaje es habitual usar una clase similar a ExpressionSet llamada DGEList” pensadas para manejar datos de contajes , definida en el paquete edgeR. Esta clase, más simple que las anteriores, utiliza listas para almacenar recuentos de “reads” e información asociada de tecnologías de secuenciación. Puede encontrarse información al respecto en la ayuda del paquete edgeR.
dgeObj <- DGEList(counts = counts.keep, lib.size = colSums(counts.keep), norm.factors = rep(1,
ncol(counts.keep)), samples = selectedTargets, group = selectedTargets$group,
genes = rownames(counts.keep), remove.zeros = FALSE)
show(dgeObj)An object of class "DGEList"
$counts
COV145 COV147 COV149 COV150 COV151
ENSG00000227232 0.06644967 0.00000000 0.2218932 0.84819243 0.07320991
ENSG00000278267 0.19934902 0.07170342 0.1479288 0.00000000 0.21962973
COV152 COV153 COV154 COV155 COV156 COV157
ENSG00000227232 0.15060933 1.3880853 0.41957493 1.163313 0.0000000 0.11736232
ENSG00000278267 0.00000000 0.0000000 0.06992916 0.174497 0.1279917 0.11736232
COV175 COV176 COV177 COV178 COV179 COV180
ENSG00000227232 0.0000000 0.9503052 1.9885554 0.1726206 0.14628246 0.6558384
ENSG00000278267 0.2236183 0.2193012 0.4349965 0.0000000 0.21942369 0.0000000
HEA061 HEA062 HEA063 HEA064 HEA065 HEA066
ENSG00000227232 0.3528731 0.0000000 0.0000000 0.06183451 0.10298192 0.05475139
ENSG00000278267 0.2352487 0.2880235 0.1926772 0.49467608 0.15447288 0.49276252
HEA067 HEA068 HEA069 HEA070 HEA071 HEA181
ENSG00000227232 0.0000000 0.06713275 0.06693886 0.04887949 0.1122282 0.2003508
ENSG00000278267 0.4268581 0.13426549 0.46857204 0.09775898 0.3366846 0.2003508
HEA182 HEA183 HEA184 HEA185 HEA186
ENSG00000227232 0.28009698 0.05465108 0.0000000 0.0000000 0.41517606
ENSG00000278267 0.14004849 0.32790646 0.6985910 0.3036215 0.17793260
[ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...
$samples
group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1 COV145 COVID red
COV147 COVID 999530.1 1 COV147 COVID red
COV149 COVID 999667.8 1 COV149 COVID red
COV150 COVID 999872.5 1 COV150 COVID red
COV151 COVID 999904.2 1 COV151 COVID red
29 more rows ...
$genes
genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
33034 more rows ...Aunque podríamos haber creado el objeto a partir de todas las muestras, y haber realizado la extracción de genes y muestras posteriormente, hemos optado por no hacerlo para facilitar el seguimiento del proceso.
3.4 Normalización
Además de estandarizar los contajes, es importante eliminar otros sesgos de composición entre librerías. Esto puede hacerse aplicando la normalización por el método TMM que genera un conjunto de factores de normalización, tal que producto de estos factores y los tamaños de librería (el número de secuencias de cada muestra) definen el tamaño efectivo de dichas muestras, es decir el peso real que se les asignará en las comparaciones posteriores.
Aunque esto puede parecer artificial, no lo es porque la normalización tiene en cuenta otros factores, como el sesgo de composición entre librerías, que podrían hacer que los mismos valores en distintas muestras no reflejaran su importancia relativa.
La función calcNormFactors, de la librería edgeR, calcula los factores de normalización mencionados.
library(edgeR)
dgeObj_norm <- calcNormFactors(dgeObj)
head(dgeObj_norm$samples, 10) group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1.0018278 COV145 COVID red
COV147 COVID 999530.1 1.1010464 COV147 COVID red
COV149 COVID 999667.8 0.9553476 COV149 COVID red
COV150 COVID 999872.5 0.9479175 COV150 COVID red
COV151 COVID 999904.2 0.8996072 COV151 COVID red
COV152 COVID 999763.7 1.0033710 COV152 COVID red
COV153 COVID 999874.5 0.9870834 COV153 COVID red
COV154 COVID 999736.9 0.8538699 COV154 COVID red
COV155 COVID 999399.0 0.8230210 COV155 COVID red
COV156 COVID 999868.7 0.8878678 COV156 COVID redEsto no modificará la matriz de contajes, pero actualizará los factores de normalización en el objeto DGEList (sus valores predeterminados son 1).
head(dgeObj$samples) group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1 COV145 COVID red
COV147 COVID 999530.1 1 COV147 COVID red
COV149 COVID 999667.8 1 COV149 COVID red
COV150 COVID 999872.5 1 COV150 COVID red
COV151 COVID 999904.2 1 COV151 COVID red
COV152 COVID 999763.7 1 COV152 COVID redhead(dgeObj_norm$samples) group lib.size norm.factors samples group.1 cols
COV145 COVID 999825.8 1.0018278 COV145 COVID red
COV147 COVID 999530.1 1.1010464 COV147 COVID red
COV149 COVID 999667.8 0.9553476 COV149 COVID red
COV150 COVID 999872.5 0.9479175 COV150 COVID red
COV151 COVID 999904.2 0.8996072 COV151 COVID red
COV152 COVID 999763.7 1.0033710 COV152 COVID redEs decir, aunque no se observen cambios en la matriz de contajes, cuando se utilizan estos factores de normalización en algún cálculo la importancia de las distintas columnas se tendrá en cuenta.
3.5 Transformación logarítmica
Las transformaciones anteriores buscan compensar el tamaño distinto de las librerías o la distinta composición de éstas, pero las distribuciones de los contajes en cada muestra son asimétricas.
boxplot(dgeObj_norm$counts, col = dgeObj_norm$samples$cols, las = 2, cex.axis = 0.7,
main = "Contajes normalizados", ylim = c(0, 10000))
Para finalizar el preprocesado se toman logaritmo base dos de los contajes
log2count_norm <- cpm(dgeObj_norm, log = TRUE)boxplot(log2count_norm, col = dgeObj_norm$samples$cols, las = 2, cex.axis = 0.7,
main = "Contajes normalizados (log2))")
source("./Rcode/niceBoxPlot.R")
niceBoxPlot(log2count_norm)
Esta será nuestra matriz de partida para los análisis siguientes
4 Exploración de los datos
Una vez descartados los genes poco expresados y con los recuentos almacenados en un objeto DGEList, podemos`proceder a realizar algunos gráficos exploratorios para determinar si los datos aparentan buena calidad y/o si presentan algun problema.
4.1 Distribución de los contajes
Un diagrama de cajas con los datos, normalizados o no, muestra que la distribución de los contajes es muy asimétrica, lo que justifica la decisión de trabajar con los logaritmos de los datos.
La transformación logarítmica puede hacerse directamente pero es mejor usar la función cpm, como se ha hecho, que agrega una pequeña cantidad para evitar tomar logaritmos de cero.
par(mfrow = c(2, 1))
rawCounts <- dgeObj_norm$counts
boxplot(rawCounts, ylab = "CPM", las = 2, xlab = "", col = dgeObj$samples$cols, cex.axis = 0.6,
main = "Distribución de contajes")
boxplot(log2count_norm, ylab = "Log2-CPM", las = 2, xlab = "", col = dgeObj$samples$cols,
cex.axis = 0.6, main = "Distribución de log(contajes)")
abline(h = median(log2count_norm), col = "blue")
par(mfrow = c(1, 1))4.2 Análisis de similaridad entre las muestras
4.2.1 Distancia entre muestras
La función dist permite calcular una matriz de distancias que contiene las comparaciones dos a dos entre todas las muestras. Por defecto se utiliza una distancia euclídea.
sampleDists <- dist(t(log2count_norm))
round(sampleDists, 1) COV145 COV147 COV149 COV150 COV151 COV152 COV153 COV154 COV155 COV156
COV147 80.0
COV149 70.7 82.0
COV157 COV175 COV176 COV177 COV178 COV179 COV180 HEA061 HEA062 HEA063
COV147
COV149
HEA064 HEA065 HEA066 HEA067 HEA068 HEA069 HEA070 HEA071 HEA181 HEA182
COV147
COV149
HEA183 HEA184 HEA185
COV147
COV149
[ reached getOption("max.print") -- omitted 31 rows ]Las matrices de distancias se pueden visualizar directamente con un heatmap.
library(factoextra)
fviz_dist(sampleDists)
Como puede verse _las muestras tienden a agruparse por el factor SANO/COVID, aunque una de las muestras. COV155 se separa del resto de las del grupo COVID.
4.2.2 Agrupamiento jerárquico
Un agrupamiento jerárquico proporciona una representación alternativa, también basada en la matriz de distancias.
hc <- hclust(sampleDists)
plot(hc, labels = colnames(log2count_norm), main = "Agrpamiento jerárquico de las muestras",
cex = 0.8)
source("./Rcode/niceDendrogram.R")
hc <- hclust(sampleDists, method = "ward.D2")
niceDendrogram(hc)
Una de las muestras COVID parece más similar a las saludables que a las otras COVID.
4.2.3 Visualización en dimensión reducida
col.status <- dgeObj_norm$samples$cols
plotMDS(log2count_norm, col = col.status, main = "Status", cex = 0.7)
source("./Rcode/niceMDS.R")
niceMDS(dgeObj_norm)
Como puede verse, el gráfico muestra la misma agrupación “natural” y el mismo comportamiento atípico de una muestra.
5 Análisis de expresión diferencial
El objetivo del análisis de expresión diferencial es seleccionar genes cuya expresión difiere entre grupos.
5.1 Selección de genes usando limma-Voom
La ventaja principal de esta aproximación es que permite trabajar con toda la flexibilidad de los modelos lineales para representar diseños experimentales, y, en muchos casos , aprovechar la experiencia previa del usuario en el manejo de limma.
5.1.1 Matriz de diseño y de contrastes
Utilizando la variable group podemos definir una matriz de diseño y, sobre ésta, los contrastes que nos interesan.
group = as.factor(dgeObj_norm$samples$group)
design = model.matrix(~0 + group)
colnames(design) = gsub("group", "", colnames(design))
row.names(design) = sampleNames
design COVID SANO
COV145 1 0
COV147 1 0
COV149 1 0
COV150 1 0
COV151 1 0
COV152 1 0
COV153 1 0
COV154 1 0
COV155 1 0
COV156 1 0
COV157 1 0
COV175 1 0
COV176 1 0
COV177 1 0
COV178 1 0
COV179 1 0
COV180 1 0
HEA061 0 1
HEA062 0 1
HEA063 0 1
HEA064 0 1
HEA065 0 1
HEA066 0 1
HEA067 0 1
HEA068 0 1
HEA069 0 1
HEA070 0 1
HEA071 0 1
HEA181 0 1
HEA182 0 1
HEA183 0 1
HEA184 0 1
HEA185 0 1
HEA186 0 1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"Dado que estamos interesados en las diferencias entre los grupos, necesitamos especificar qué comparaciones queremos llevar a cabo. Las comparaciones de interés se puede especificar utilizando la función makeContrasts. La matriz de contraste indica qué columnas de la matriz design vamos a comparar. En este caso tan sólo se llevará a cabo una comparación.
cont.matrix = makeContrasts(CONTROLvsCOVID = COVID - SANO, levels = colnames(design))
cont.matrix Contrasts
Levels CONTROLvsCOVID
COVID 1
SANO -15.1.2 Transformación de los contajes
voomObj <- voom(dgeObj_norm, design)
voomObjAn object of class "EList"
$genes
genes
ENSG00000227232 ENSG00000227232
ENSG00000278267 ENSG00000278267
ENSG00000233750 ENSG00000233750
ENSG00000268903 ENSG00000268903
ENSG00000269981 ENSG00000269981
33034 more rows ...
$targets
group lib.size norm.factors samples group.1 cols
COV145 COVID 1001653.3 1.0018278 COV145 COVID red
COV147 COVID 1100529.0 1.1010464 COV147 COVID red
COV149 COVID 955030.2 0.9553476 COV149 COVID red
COV150 COVID 947796.6 0.9479175 COV150 COVID red
COV151 COVID 899521.0 0.8996072 COV151 COVID red
29 more rows ...
$E
COV145 COV147 COV149 COV150 COV151
ENSG00000227232 -0.8223650 -1.1381984 -0.4037625 0.50837551 -0.6500950
ENSG00000278267 -0.5183002 -0.9448596 -0.5597126 -0.92265092 -0.3219038
COV152 COV153 COV154 COV155 COV156
ENSG00000227232 -0.6246522 0.9358596 0.1073289 1.0159241 -0.8282289
ENSG00000278267 -1.0045156 -0.9810644 -0.5828556 -0.2862517 -0.4994115
COV157 COV175 COV176 COV177 COV178
ENSG00000227232 -0.74055345 -0.9859431 0.5497283 1.3062240 -0.6882933
ENSG00000278267 -0.74055345 -0.4526424 -0.4619603 -0.1060517 -1.1161581
COV179 COV180 HEA061 HEA062 HEA063
ENSG00000227232 -0.654532686 0.1342634 -0.20656331 -1.01143480 -1.1012204
ENSG00000278267 -0.499855856 -1.0746764 -0.42066201 -0.35512422 -0.6309653
HEA064 HEA065 HEA066 HEA067 HEA068
ENSG00000227232 -0.9265423 -0.7259291 -0.78574256 -1.1218171 -0.9590727
ENSG00000278267 -0.1024608 -0.6077105 0.05386471 -0.2313968 -0.7976723
HEA069 HEA070 HEA071 HEA181 HEA182
ENSG00000227232 -1.0099326 -0.8215555 -0.6475041 -0.4772445 -0.4479095
ENSG00000278267 -0.2372664 -0.6984811 -0.1968898 -0.4772445 -0.7333818
HEA183 HEA184 HEA185 HEA186
ENSG00000227232 -0.8517698 -1.0495769 -0.9932084 -0.1628122
ENSG00000278267 -0.2738825 0.2117625 -0.3086203 -0.5957197
[ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...
$weights
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 6.976324 6.423290 7.345017 7.414329 7.913732 6.966953 7.070951
[2,] 13.344117 9.792171 16.233083 16.785098 21.291716 13.276959 14.091120
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 8.568175 9.129126 8.046130 6.792973 7.049157 7.046103 6.946915
[2,] 27.856177 34.175509 22.694541 12.084889 13.891449 13.863695 13.134382
[,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,] 6.497274 6.878669 6.667025 21.376991 19.139568 14.546448 14.827401
[2,] 10.280091 12.659254 11.283734 8.528972 8.192672 7.560903 7.603134
[,22] [,23] [,24] [,25] [,26] [,27] [,28]
[1,] 20.095913 24.685105 13.715056 13.110198 11.661751 22.943084 24.332783
[2,] 8.340050 8.955089 7.428218 7.308304 7.031769 8.740787 8.912791
[,29] [,30] [,31] [,32] [,33] [,34]
[1,] 22.360747 15.055081 19.758022 16.994284 20.284112 17.780471
[2,] 8.666115 7.636840 8.288627 7.906856 8.368394 8.008534
[ reached getOption("max.print") -- omitted 3 rows ]
33034 more rows ...
$design
COVID SANO
COV145 1 0
COV147 1 0
COV149 1 0
COV150 1 0
COV151 1 0
29 more rows ...5.1.3 Selección de genes diferencialmente expresados
Como en el caso de los microarrays el objeto voomObj y las matrices de diseño y contrastes se utilizaran para ajustar un modelo y, a continuación realizar las comparaciones especificadas sobre el modelo ajustado. El proceso finaliza con la regularización del estimador del error usando la función eBayes.
fit <- lmFit(voomObj)
fit.cont <- contrasts.fit(fit, cont.matrix)
fit.cont <- eBayes(fit.cont)5.1.4 Top tables
Los resultados de un análisis de expresión diferencial se pueden extraer con la función topTable. Esta función genera una tabla de resultados cuyas columnas contienen información acerca de los genes y la diferencia entre los grupos comparados. Concretamente:
toptab <- topTable(fit.cont, coef = 1, sort.by = "p", number = nrow(fit.cont))
head(toptab) genes logFC AveExpr t P.Value
ENSG00000167900 ENSG00000167900 2.832423 4.111065 12.36669 8.088344e-15
ENSG00000090889 ENSG00000090889 2.476447 1.292066 11.83766 3.015984e-14
ENSG00000111206 ENSG00000111206 2.423923 2.018435 11.80083 3.309678e-14
ENSG00000166851 ENSG00000166851 3.223189 2.615540 11.65570 4.781203e-14
ENSG00000135476 ENSG00000135476 2.236975 1.705254 11.53985 6.425179e-14
ENSG00000123485 ENSG00000123485 2.553959 1.914888 11.46276 7.828950e-14
adj.P.Val B
ENSG00000167900 2.672308e-10 23.54910
ENSG00000090889 3.644948e-10 22.09352
ENSG00000111206 3.644948e-10 21.99541
ENSG00000166851 3.949154e-10 21.70990
ENSG00000135476 4.013294e-10 21.35292
ENSG00000123485 4.013294e-10 21.189535.1.5 Visualización de los resultados
5.1.5.1 Significación bioogica vs Significación Estadística
Para visualizar los resultados podemos usar un volcanoPlot:
volcanoplot(fit.cont, coef = 1, highlight = 100, main = "COVID vs SANO")
Es posible también crear Volcanoplots interactivos
source("./Rcode/niceInteractiveVolcano.R")
niceInteractiveVolcano(fit.cont)5.1.5.2 Perfiles de expresión
Con el fin de observar si existen perfiles de expresión diferenciados podemo realizar un mapa de colores con los genes más diferencialmente expresados.
Es decir, fijamos un criterio de selección de genes y retenemos aquellos componentes de la tabla de resultados que lo cumplen. Por ejemplo: Genes con un p-balor ajustado inferior a 0.001 y un `fold-change’ superior a 2.
topGenesBas <- rownames(subset(toptab, (abs(logFC) > 2) & (adj.P.Val < 0.01)))
length(topGenesBas)[1] 199Con la matriz de expresión de los genes que verifican dicha condición se puede construir un heatmap.
library(pheatmap)
mat <- log2count_norm[topGenesBas, ]
mat <- mat - rowMeans(mat)
pheatmap(mat)
También es posible crear un Heatmap interactivo:
source("./Rcode/niceInteractiveHeatMap.R")
mat <- log2count_norm[topGenesBas, ]
mat <- mat - rowMeans(mat)
niceInteractiveHeatMap(mat)Los dos grupos estan diferenciados, sobre todo en un subconjunto de genes en donde las expresiones toman signos (colores) distintos. En otro de los grupos parece que, en el grupo de individuos sanos los genes diferencialmente expresados se encuentran sobre-expresados en el grupo COVID, y apenas expresados en el grupo sanos.
5.2 Análisis de expresión diferencial usando el paquete edgeR
El análisis con edgeR es similar al anterior (se originan en el mismo equipo de investigación) pero la modelización es distinta.
El análisis utiliza un GLM pero, en una forma que recuerda a lo que se hace con limma, realiza un paso adicional,en el que se calcula una estimación mejorada de la dispersión (variabilidad) de las muestras que integra las estimaciones individuales y la global mediante estimación Bayes empírica.
y = estimateDisp(dgeObj_norm, design, robust = TRUE)
plotBCV(y)
Con este objeto, que añade a los contajes normalizados los estimadores mejorados de dispersión, se ajusta un modelo lineal generalizado con distribución binomial para los errores.
fit <- glmQLFit(y, design, robust = TRUE)Una vez ajustado el modelo se procede a construir el contraste y realizar el test.
De hecho podemos usar la matriz de contrastes que construímos para limma voom, en la misma forma que hemos reutilizado la de diseño.
res <- glmQLFTest(fit, contrast = cont.matrix)Los resultados se almacenan en un objeto similar a la ’ topTable’ de limma.
topTags_edge <- topTags(res, n = dim(log2count_norm)[1]) # todos los genes
head(topTags_edge)Coefficient: 1*COVID -1*SANO
genes logFC logCPM F PValue
ENSG00000123485 ENSG00000123485 2.969009 2.874508 143.3978 6.410064e-15
ENSG00000167900 ENSG00000167900 3.052165 4.946915 165.9368 1.024164e-14
ENSG00000187837 ENSG00000187837 1.719557 4.813853 135.1654 1.889383e-14
ENSG00000111206 ENSG00000111206 2.834520 2.919472 132.0747 2.767000e-14
ENSG00000178999 ENSG00000178999 2.987952 3.225967 135.0640 4.056862e-14
ENSG00000166851 ENSG00000166851 3.638997 3.745942 144.7154 5.207210e-14
FDR
ENSG00000123485 1.691868e-10
ENSG00000167900 1.691868e-10
ENSG00000187837 2.080777e-10
ENSG00000111206 2.285473e-10
ENSG00000178999 2.680693e-10
ENSG00000166851 2.867350e-10Podemos seleccionar los genes más diferencialmente expresados de la misma forma que hicimos con limma-voom
topGenes_edge <- rownames(subset(topTags_edge$table, (abs(logFC) > 2) & (FDR < 0.01)))
length(topGenes_edge)[1] 376Obsérvese que la lista de genes seleccionados es muy similar en ambos casos.
library(ggvenn)
x = list(LimmaVoom = topGenesBas, edgeR = topGenes_edge)
ggvenn(x, fill_color = c("#0073C2FF", "#EFC000FF"), stroke_size = 0.5, set_name_size = 3)
6 Anotación de resultados y análisis de significación biológica
Para el análisis de significación se utilizan dos listas de transcritos:
- La lista de transcritos diferencialmenete expresados
- La lista de todos los tránscritos o “Universo”
topGenes <- union(topGenesBas, topGenes_edge)
length(topGenes)[1] 383universe <- rownames(toptab)
length(universe)[1] 330396.1 Anotación de los identificadores
Esto es posible, y de hecho sencillo de llevar a cabo, usando el paquete annotate.
library(org.Hs.eg.db)
AnnotationDbi::keytypes(org.Hs.eg.db) [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
[6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[11] "GENETYPE" "GO" "GOALL" "IPI" "MAP"
[16] "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM"
[21] "PMID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
[26] "UNIPROT" topAnots = AnnotationDbi::select(org.Hs.eg.db, topGenes, c("SYMBOL", "ENTREZID",
"GENENAME"), keytype = "ENSEMBL")
head(topAnots) ENSEMBL SYMBOL ENTREZID GENENAME
1 ENSG00000167900 TK1 7083 thymidine kinase 1
2 ENSG00000090889 KIF4A 24137 kinesin family member 4A
3 ENSG00000111206 FOXM1 2305 forkhead box M1
4 ENSG00000166851 PLK1 5347 polo like kinase 1
5 ENSG00000135476 ESPL1 9700 extra spindle pole bodies like 1, separase
6 ENSG00000123485 HJURP 55355 Holliday junction recognition proteindim(topAnots)[1] 385 4Como puede verse, el número de anotaciones es el mismo que el de identificadores ENSEMBL, lo que podría llevar a pensar que, es posible que, antes de subir los datos a GEO se hayan agrupado los contajes por genes.
Para la anotación del universo se procederá igual.
univAnots = AnnotationDbi::select(org.Hs.eg.db, universe, c("SYMBOL", "ENTREZID",
"GENENAME"), keytype = "ENSEMBL")
head(univAnots) ENSEMBL SYMBOL ENTREZID GENENAME
1 ENSG00000167900 TK1 7083 thymidine kinase 1
2 ENSG00000090889 KIF4A 24137 kinesin family member 4A
3 ENSG00000111206 FOXM1 2305 forkhead box M1
4 ENSG00000166851 PLK1 5347 polo like kinase 1
5 ENSG00000135476 ESPL1 9700 extra spindle pole bodies like 1, separase
6 ENSG00000123485 HJURP 55355 Holliday junction recognition proteindim(univAnots)[1] 33202 4En este caso se observa como hay más anotaciones que transcritos, lo que sugiere que múltiples transcritos han sido mapeados en el mismo gen.
6.2 Análisis de enriquecimiento
El paquete clusterProfiler admite identificadores de tipo ENSEMBL y permite gran variedad de análisis complementarios al enriquecimiento por lo que, es una de las mejores opciones para el análisis de significación biológica.
library(clusterProfiler)
library(org.Hs.eg.db)
ego = enrichGO(gene = topGenes, universe = universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db,
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable = TRUE)
head(ego[, -8], n = 5) ID
GO:0006958 GO:0006958
GO:0002455 GO:0002455
GO:0006956 GO:0006956
GO:0016064 GO:0016064
GO:0019724 GO:0019724
Description
GO:0006958 complement activation, classical pathway
GO:0002455 humoral immune response mediated by circulating immunoglobulin
GO:0006956 complement activation
GO:0016064 immunoglobulin mediated immune response
GO:0019724 B cell mediated immunity
GeneRatio BgRatio pvalue p.adjust qvalue Count
GO:0006958 80/342 120/16130 4.046401e-107 1.215134e-103 1.117659e-103 80
GO:0002455 81/342 135/16130 9.256747e-103 1.389901e-99 1.278405e-99 81
GO:0006956 81/342 148/16130 4.506649e-98 4.511156e-95 4.149280e-95 81
GO:0016064 81/342 200/16130 3.900362e-84 2.928197e-81 2.693303e-81 81
GO:0019724 81/342 203/16130 1.728105e-83 1.037900e-80 9.546418e-81 81Con los resultados del análisis de enriquecimiento se pueden llevar a cabo distintas visualizaciones cuya interpretación exacta puede verse en el manual de clusterProfiler
dotplot(ego, showCategory = 7)
library(ggplot2)
ego2 = simplify(ego)
cnetplot(ego2, showCategory = 3, cex_category = 0.3, cex_label_category = 0.7, cex_gene = 0.2,
cex_label_gene = 0.4, circular = TRUE, colorEdge = TRUE)
library(enrichplot)
goplot(ego2, showCategory = 6, cex = 0.1)
heatplot(ego2)
term_similarity_matrix = pairwise_termsim(ego)
emapplot(term_similarity_matrix, showCategory = 15, group_category = TRUE, group_legend = TRUE)